home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus Leser 15 / Amiga Plus Leser CD 15.iso / Tools / Development / yacas_alg / yacas_morphos / share / yacas / include / anumber.inl < prev    next >
Encoding:
Text File  |  2002-03-13  |  15.4 KB  |  619 lines
  1.  
  2.  
  3. /*
  4. #include <stdio.h>
  5. template<class T>
  6. inline void Thunk(T& a)
  7. {
  8.   int i;
  9.   for (i=0;i<a.NrItems();i++)
  10.     printf("%d ",a[i]);
  11.   printf("\n");  
  12. }
  13. */
  14.  
  15. /* BaseTimesInt : multiply a with one digit in the range 0..(aBase-1)
  16.  */
  17. template<class T>
  18. inline void BaseTimesInt(T& a,PlatDoubleWord aNumber, PlatDoubleWord aBase)
  19. {
  20. //    if (a[a.NrItems()-1] != 0)
  21. //      a.Append(0);
  22.  
  23.     PlatDoubleWord carry=0;
  24.     LispInt i;
  25.     LispInt nr=a.NrItems();
  26.  
  27.     typename T::ElementTypePtr aptr = &a[0];
  28.     for (i=0;i<nr;i++)
  29.     {
  30.         PlatDoubleWord word = ((PlatDoubleWord)(*aptr))*aNumber+carry;
  31.         PlatWord digit = (PlatWord)(word % aBase);
  32.         PlatWord newCarry= (PlatWord)(word / aBase);
  33.         *aptr = digit;
  34.         aptr++;
  35.         carry= newCarry;
  36.     }
  37.     if (carry)
  38.     {
  39.       a.Append((typename T::ElementType)carry);
  40.       carry = 0;
  41.     }
  42.     LISPASSERT(carry == 0);
  43. }
  44.  
  45. template<class T>
  46. inline void BaseDivideInt(T& a,PlatDoubleWord aNumber, PlatDoubleWord aBase, PlatDoubleWord& aCarry)
  47. {
  48. //    if (a[a.NrItems()-1] != 0)
  49.  
  50.     PlatDoubleWord carry=0;
  51.     LispInt i;
  52.     LispInt nr=a.NrItems();
  53.  
  54.     typename T::ElementTypePtr aptr = &a[0];
  55.     for (i=nr-1;i>=0;i--)
  56.     {
  57.         PlatDoubleWord word = ((carry*aBase)+((PlatDoubleWord)(aptr[i])));
  58.         PlatWord digit = (PlatWord)(word / aNumber);
  59.         PlatWord newCarry= (PlatWord)(word % aNumber);
  60.         aptr[i] = digit;
  61.         carry= newCarry;
  62.     }
  63.     //carry now is the remainder
  64.     aCarry = carry;
  65. }
  66.  
  67.  
  68. /* GrowDigits : add digits to a until it has aDigits digits
  69.  */
  70. template<class T>
  71. inline void GrowDigits(T& a,LispInt aDigits)
  72. {
  73.     LispInt i;
  74.  
  75.     if (aDigits <= a.NrItems())
  76.         return;
  77.     
  78.     /*
  79.      LispInt nrToAdd = aDigits-a.NrItems();
  80.  
  81.     for (i=0;i<nrToAdd;i++)
  82.         a.Append(0);
  83. */
  84.     LispInt origSize = a.NrItems();
  85.     a.GrowTo(aDigits);
  86.     a.SetNrItems(aDigits);
  87.     if (aDigits<=origSize)
  88.         return;
  89.     typename T::ElementType* ptr = &a[origSize];
  90.     for (i=origSize;i<aDigits;i++)
  91.         *ptr++ = 0;
  92. }
  93.  
  94. /* BaseAdd : destructively add aSource to aTarget, in base aBase.
  95.  */
  96. template<class T>
  97. inline void BaseAdd(T& aTarget, const T& aSource, LispInt aBase)
  98. {
  99.     // Initialize result
  100.  
  101.     GrowDigits(aTarget,aSource.NrItems());
  102.     aTarget.Append(0);
  103.     
  104.     LispInt nr1 = aTarget.NrItems();
  105.     LispInt nr2 = aSource.NrItems();
  106.     LispInt nr;
  107.  
  108.     // nr represents min(nr1,nr2), the number of digits to add
  109.     if (nr1>nr2)
  110.         nr=nr2;
  111.     else
  112.         nr=nr1;
  113.     
  114.     PlatDoubleWord carry=0;
  115.     LispInt digit;
  116.  
  117.    typename T::ElementTypePtr sourcePtr = &aSource[0];
  118.    typename T::ElementTypePtr targetPtr = &aTarget[0];
  119.    for (digit=0;digit<nr;digit++)
  120.     {
  121.         PlatDoubleWord word;
  122.         word = (PlatDoubleWord)targetPtr[digit] +
  123.             (PlatDoubleWord)sourcePtr[digit] + carry;
  124.          PlatDoubleWord newDigit = (word%aBase);
  125.          PlatDoubleWord newCarry = (word/aBase);
  126.          targetPtr[digit] = (typename T::ElementType)newDigit;
  127.          carry          = newCarry;
  128.     }
  129.     while (carry != 0)
  130.     {
  131.         PlatSignedDoubleWord ww = targetPtr[nr];
  132.         ww+=carry;
  133.         targetPtr[nr] = ww%aBase;
  134.         carry = ww/aBase;
  135.         nr++;
  136.     }
  137. }
  138.  
  139. template<class T>
  140. inline void BaseSubtract(T& aResult, T& a2, LispInt offset)
  141. {
  142.     if (IsZero(a2))
  143.         return;
  144.     LISPASSERT(!IsZero(a2));
  145.     // Initialize result
  146.     LispInt nr = a2.NrItems();
  147.  
  148.     typename T::ElementTypePtr resultPtr = &aResult[0];
  149.     typename T::ElementTypePtr a2ptr = &a2[0];
  150.  
  151.     while (a2ptr[nr-1] == 0)
  152.         nr--;
  153.  
  154.     // Subtract on a per-digit basis
  155.     PlatSignedDoubleWord carry=0;
  156.     LispInt digit;
  157.  
  158.     for (digit=0;digit<nr;digit++)
  159.     {
  160.         PlatSignedDoubleWord word;
  161.         word = ((PlatSignedDoubleWord)resultPtr[digit+offset]) -
  162.             ((PlatSignedDoubleWord)a2ptr[digit]) +
  163.             (PlatSignedDoubleWord)carry;
  164.         carry=0;
  165.         while (word<0)
  166.         {
  167.             word+=WordBase;
  168.             carry--;
  169.         }
  170.         resultPtr[digit+offset] = ((PlatWord)(word%WordBase));
  171.     }
  172.  
  173.     while (carry != 0)
  174.     {
  175.         LISPASSERT(nr+offset<aResult.NrItems());
  176.  
  177.         LispInt newCarry = 0;
  178.         PlatSignedDoubleWord ww = resultPtr[nr+offset]+carry;
  179.         while (ww<0)
  180.         {
  181.             ww = ww + WordBase;
  182.             newCarry = newCarry - 1;
  183.         }
  184.         resultPtr[nr+offset]=(typename T::ElementType)ww;
  185.         carry = newCarry;
  186.         offset++;
  187.     }
  188. }
  189.  
  190. /* BaseIntNumber : convert a number into a different base,
  191.  * using growing arrays.
  192.  */
  193. template<class T>
  194. inline void BaseIntNumber(T& aTarget, LispInt aNumber, PlatWord aBase)
  195. {
  196.     aTarget.SetNrItems(0);
  197.     while (aNumber != 0)
  198.     {
  199.         LispInt digit = aNumber%aBase;
  200.         aTarget.Append(digit);
  201.         aNumber/=aBase;
  202.     }
  203.     if (aTarget.NrItems() == 0)
  204.         aTarget.Append(0);
  205. }
  206.  
  207. // BaseAddMultiply : multiply x and y, and add result to aTarget
  208. //
  209. template<class T>
  210. inline void BaseAddMultiply(T& aTarget, T& x, T& y, LispInt aBase)
  211. {
  212.     LispInt nrx=x.NrItems();
  213.     LispInt nry=y.NrItems();
  214.     GrowDigits(aTarget,nrx+nry+1);
  215.     LispInt ix,iy;
  216.  
  217.     typename T::ElementType *targetPtr = &aTarget[0];
  218.     typename T::ElementType *xPtr = &x[0];
  219.     typename T::ElementType *yPtr = &y[0];
  220.     for (ix=0;ix<nrx;ix++)
  221.     {
  222.         PlatDoubleWord carry = 0;
  223.         for (iy=0;iy<nry;iy++)
  224.         {
  225.             PlatDoubleWord word =
  226.                 static_cast<PlatDoubleWord>(targetPtr[ix+iy])+
  227.                 static_cast<PlatDoubleWord>(xPtr[ix])*
  228.                 static_cast<PlatDoubleWord>(yPtr[iy])+carry;
  229.             // This calculates aTarget[ix+iy]+x[ix]*y[iy]+carry;
  230.  
  231.  
  232.             targetPtr[ix+iy] = (typename T::ElementType)(word % aBase);
  233.             carry          = word / aBase;
  234.         }
  235.         targetPtr[ix+nry] += (typename T::ElementType)(carry);
  236.     }
  237. }
  238.  
  239.  
  240. #ifdef USE_KARATSUBA
  241.  
  242. const int cKaratsubaCutoff = 8;
  243.  
  244. /*
  245. #include <stdio.h>
  246.  
  247. inline void BaseKaratsubaMultiply(ANumber& aTarget, ANumber& x, ANumber& y, LispInt aSize, LispInt aBase)
  248. {
  249.     LispInt halfSize = aSize/2;
  250.  
  251.     ANumber F0;
  252.     ANumber F1;
  253.     ANumber G0;
  254.     ANumber G1;
  255.  
  256.     GrowDigits(F0,halfSize);
  257.     GrowDigits(F1,halfSize);
  258.     GrowDigits(G0,halfSize);
  259.     GrowDigits(G1,halfSize);
  260.     LispInt i;
  261.     for (i=0;i<halfSize;i++)
  262.     {
  263.         F0[i] = x[i];
  264.         F1[i] = x[halfSize+i];
  265.         G0[i] = y[i];
  266.         G1[i] = y[halfSize+i];
  267.     }
  268.     
  269.     ANumber F0G0;
  270.     ANumber F1G1;
  271.     ANumber Fsum,Gsum;
  272.     ANumber combined;
  273.  
  274.     Multiply(F0G0, F0, G0);
  275.     Multiply(F1G1, F1, G1);
  276.     Add(Fsum,F0,F1);
  277.     Add(Gsum,G0,G1);
  278.     Multiply(combined, Fsum, Gsum);
  279.     ANumber s1;
  280.     Subtract(s1,combined,F0G0);
  281.     ANumber s2;
  282.     Subtract(s2,s1,F1G1);
  283.  
  284.     BaseShiftLeft(F1G1, aSize   *(8*sizeof(PlatWord)));
  285.     BaseShiftLeft(s2  , halfSize*(8*sizeof(PlatWord)));
  286.  
  287.     ANumber s3;
  288.     Add(s3,F1G1,s2);
  289.     Add(aTarget,s3,F0G0);
  290. }
  291.  
  292. inline void BaseAddMultiplyK(ANumber& aTarget, ANumber& x, ANumber& y, LispInt aBase)
  293. {
  294.     LispInt nrx=x.NrItems();
  295.     LispInt nry=y.NrItems();
  296.     LispInt i, i2;
  297.  
  298.     // Too small?  Use O(n^2) multiplication
  299.     if (nrx+nry+1 <= cKaratsubaCutoff)
  300.     {
  301.         BaseAddMultiply(aTarget, x, y, aBase);
  302.         return;
  303.     }
  304.  
  305.     LispInt max = (nrx > nry) ? nrx : nry;
  306.     LispInt maxtwos;
  307.     for (maxtwos = 1; maxtwos < max; maxtwos *= 2);
  308.  
  309.     GrowDigits(aTarget,maxtwos);
  310.     GrowDigits(x,maxtwos);
  311.     GrowDigits(y,maxtwos);
  312.     BaseKaratsubaMultiply(aTarget, x, y, maxtwos, aBase);
  313. }
  314. */
  315.  
  316.  
  317. /* BaseAddMultiplyK : multiply x and y, and add result to aTarget
  318.  *                      using Karatsuba multiplication.  This function
  319.  *                      just makes the numbers passed into even ones,
  320.  *                      by adding zeros to the front, then multiplies.
  321.  *
  322.  *                      Needs enough memory to hold 4*(smallest power
  323.  *                      of 2 bigger than operands) digits.
  324.  *
  325.  *                      Code based on
  326.  *
  327.  *            Here's a crude diagram of how the numbers are stored:
  328.  *            (scratch space is where the x & y sums are stored)
  329.  *
  330.  *                 ___________ ___________ _______________________
  331.  *                |           |           |  x & y scratch space  |
  332.  *                |     x     |     y     |  starts on <-- side   |
  333.  *                |___________|___________|____ and works --> ____| 
  334.  *
  335.  */
  336.  
  337.  
  338. template<class T>
  339. inline void BaseAddMultiplyK(T& aTarget, T& x, T& y, LispInt aBase)
  340. {
  341.     LispInt nrx=x.NrItems();
  342.     LispInt nry=y.NrItems();
  343.     LispInt i, i2;
  344.     typename T::ElementType *iArray, *iXArray, *iYArray, *iSumArray;
  345.  
  346.     // Too small?  Use O(n^2) multiplication
  347.     if (nrx+nry+1 <= cKaratsubaCutoff) {
  348.         BaseAddMultiply(aTarget, x, y, aBase);
  349.         return;
  350.     }
  351.  
  352.     i = (nrx > nry) ? nrx : nry;
  353.     for (i2 = 1; i2 < i; i2 *= 2);
  354.  
  355.     // Allocate a array of the smallest power of 2 larger than both
  356.     // numbers * 4 elements, clear it to zero, then copy x into the 
  357.     // into the start, y above it, and use the high half for recursion.
  358.     // ToDo: Is this portable ?!?
  359.     iArray = (typename T::ElementType *) PlatAlloc(i2*4*sizeof(T::ElementType));
  360.     PlatMemSet((LispChar *)iArray, 0, i2*4*sizeof(T::ElementType));
  361.  
  362.     // Split our array into x and y halves
  363.     iXArray = &iArray[0];
  364.     x.CopyToExternalArray(iXArray, true);
  365.  
  366.     iYArray = &iArray[i2];
  367.     y.CopyToExternalArray(iYArray, true);
  368.  
  369.     iSumArray = &iArray[i2*2];
  370.  
  371.     BaseKaratsubaMultiply(aTarget, iXArray, iYArray, iSumArray, i2, aBase);
  372.  
  373.     PlatFree((LispChar *) iArray);
  374. }
  375.  
  376.  
  377. template<class T>
  378. void BaseKaratsubaMultiply(T& aTarget, 
  379.                            typename T::ElementType* x,
  380.                            typename T::ElementType* y,
  381.                            typename T::ElementType* aScratchSpace,
  382.                            LispInt aSize,
  383.                            LispInt aBase)
  384. {
  385.     LispInt iHalfSize = aSize/2;
  386.     LispInt i;
  387.  
  388.     typename T::ElementType *xlow= &x[0];
  389.     typename T::ElementType *xsum= &aScratchSpace[0];
  390.     typename T::ElementType *xhigh=&x[iHalfSize];
  391.  
  392.     typename T::ElementType *ylow= &y[0];
  393.     typename T::ElementType *ysum= &aScratchSpace[iHalfSize];
  394.     typename T::ElementType *yhigh=&y[iHalfSize];
  395.     T p1, p2, p3;
  396.  
  397.     // Too small?  Use O(n^2) multiplication
  398.     if (aSize+1 <= cKaratsubaCutoff)
  399.     {
  400.         T iX(x, iHalfSize), iY(y, iHalfSize);
  401.         BaseAddMultiply(aTarget, iX, iY, aBase);
  402.         return;
  403.     }
  404.  
  405.     // Compute xsum & ysum and put into the low half of the
  406.     // scratch space, the high half is used for recursion
  407.  
  408. #if 1
  409.     typename T::ElementType carryx=0;
  410.     typename T::ElementType carryy=0;
  411. #endif
  412.     for (i = 0; i < iHalfSize; i++)
  413.     {
  414. #if 1
  415.         typename T::ElementType word;
  416.         typename T::ElementType newDigit;
  417.         typename T::ElementType newCarry;
  418.  
  419.         word = (typename T::ElementType)xlow[i] +
  420.                (typename T::ElementType)xhigh[i] + carryx;
  421.         newDigit= (word%aBase);
  422.         newCarry = (word/aBase);
  423.         xsum[i] = newDigit;
  424.         carryx  = newCarry;
  425.  
  426.         word = (typename T::ElementType)ylow[i] +
  427.                (typename T::ElementType)yhigh[i] + carryy;
  428.         newDigit= (word%aBase);
  429.         newCarry = (word/aBase);
  430.         ysum[i] = newDigit;
  431.         carryy  = newCarry;
  432. #endif
  433.  
  434.         //xsum[i] = xlow[i] + xhigh[i];
  435.         //ysum[i] = ylow[i] + yhigh[i];
  436.     }
  437. #if 1
  438.     xsum[iHalfSize]+=carryx;
  439.     ysum[iHalfSize]+=carryy;
  440. #endif
  441.     
  442.     // ToDo: make array bigger, so we don't have to zero the
  443.     // scratch space after the first two calls?
  444.     //hier
  445.     BaseKaratsubaMultiply(p1, xlow, ylow, &aScratchSpace[aSize],
  446.                           iHalfSize, aBase);
  447.     PlatMemSet((LispChar *)&aScratchSpace[aSize], 0, aSize*sizeof(typename T::ElementType));
  448.     BaseKaratsubaMultiply(p2, xsum, ysum, &aScratchSpace[aSize],
  449.                           iHalfSize, aBase);
  450.     PlatMemSet((LispChar *)&aScratchSpace[aSize], 0, aSize*sizeof(typename T::ElementType));
  451.     BaseKaratsubaMultiply(p3, xhigh, yhigh, &aScratchSpace[aSize],
  452.                           iHalfSize, aBase);
  453.  
  454.     // Have a feeling these functions are slowing down the
  455.     // multipication
  456.  
  457.     BaseSubtract(p2, p1, 0);
  458.     BaseSubtract(p2, p3, 0);
  459.     BaseAdd(aTarget, p1, aBase);
  460.     BaseAdd(aTarget, p3, aBase);
  461.     BaseAdd(aTarget, p2, aBase);
  462. }
  463.  
  464. #endif
  465.  
  466. /* BaseMultiply : multiply x and y, and put result in aTarget
  467.  */
  468.  
  469. /*
  470.  #ifdef USE_KARATSUBA
  471. inline void BaseMultiply(ANumber& aTarget, ANumber& x, ANumber& y, LispInt aBase)
  472. {
  473.     aTarget.SetNrItems(1);
  474.     aTarget[0] = 0;
  475.     BaseAddMultiplyK(aTarget, x, y, aBase);
  476. }
  477. #endif
  478. */
  479.  
  480. template<class T>
  481. inline void BaseMultiply(T& aTarget, T& x, T& y, LispInt aBase)
  482. {
  483.     aTarget.SetNrItems(1);
  484.     aTarget[0] = 0;
  485. #ifdef USE_KARATSUBA
  486.     BaseAddMultiplyK(aTarget, x, y, aBase);
  487. #else
  488.     BaseAddMultiply(aTarget, x, y, aBase);
  489. #endif
  490. }
  491.  
  492.  
  493.  
  494. template<class T>
  495. inline LispBoolean IsZero(T& a)
  496. {
  497.     LispInt i;
  498.     for (i=0;i<a.NrItems();i++)
  499.         if (a[i] != 0)
  500.             return LispFalse;
  501.     return LispTrue;
  502. }
  503.  
  504.  
  505.  
  506. template<class T>
  507. inline void BaseDivide(T& aQuotient, T& aRemainder, T& a1, T& a2, PlatDoubleWord base)
  508. {
  509.     // Find the values n and m as described in Knuth II:
  510.     LispInt n,m;
  511.     n=a2.NrItems();
  512.     LISPASSERT(n>0);
  513.     LISPASSERT(a2[n-1] != 0);
  514.     
  515.     //a1.NrItems() = m+n => m = a1.NrItems()-n
  516.     m = a1.NrItems()-n;
  517.     LISPASSERT(m>=0);
  518.  
  519.     aQuotient.GrowTo(m+1);
  520.     
  521.     //D1:
  522.     //this calculates d = base/(a2[n-1]+1);
  523.     PlatDoubleWord d = base/(static_cast<PlatDoubleWord>(a2[n-1])+1);
  524.  
  525.  
  526.     BaseTimesInt(a1, d, base);
  527.     BaseTimesInt(a2, d, base);
  528.     a1.Append(0);
  529.     a2.Append(0);
  530.     
  531.     //D2:
  532.     LispInt j = m;
  533.  
  534.     while (j>=0)
  535.     {
  536.         //D3:
  537.         PlatDoubleWord q = (a1[j+n]*base+a1[j+n-1])/a2[n-1];
  538.         PlatDoubleWord r = (a1[j+n]*base+a1[j+n-1])%a2[n-1];
  539.  
  540.     REDO:
  541.         if (q == base || q*a2[n-2] > base*r+a1[j+n-2])
  542.         {
  543.             q = q - 1;
  544.             r = r + a2[n-1];
  545.             if (r < base)
  546.                 goto REDO;
  547.         }
  548.  
  549.         //D4:
  550.         ANumber sub(Precision(aQuotient));
  551.         sub.CopyFrom(a2);
  552.         BaseTimesInt(sub, q, base);
  553.         sub.Append(0);
  554.  
  555.         PlatSignedDoubleWord carry;
  556.         LispInt digit;
  557.         {//Subtract the two
  558.             //TODO this can be generalized!!!!
  559.             //
  560.             // Beware though: this is not a normal subtraction. Only a
  561.             // certain set of digits ends up being subtracted.
  562.  
  563.             // First check if qv isn't too big...
  564.             carry = 0;
  565.             for (digit=0;digit<=n;digit++)
  566.             {
  567.                 PlatSignedDoubleWord word;
  568.                 word = ((PlatSignedDoubleWord)a1[digit+j]) -
  569.                     ((PlatSignedDoubleWord)sub[digit]) +
  570.                     (PlatSignedDoubleWord)carry;
  571.                 carry=0;
  572.                 while (word<0)
  573.                 {
  574.                     word+=base;
  575.                     carry--;
  576.                 }
  577.             }
  578.             if (carry)
  579.             {
  580.                 q--;
  581.                 sub.CopyFrom(a2);
  582.                 BaseTimesInt(sub, q, base);
  583.                 sub.Append(0);
  584.             }
  585.             
  586.             carry = 0;
  587.             for (digit=0;digit<=n;digit++)
  588.             {
  589.                 PlatSignedDoubleWord word;
  590.                 word = ((PlatSignedDoubleWord)a1[digit+j]) -
  591.                     ((PlatSignedDoubleWord)sub[digit]) +
  592.                     (PlatSignedDoubleWord)carry;
  593.                 carry=0;
  594.                 while (word<0)
  595.                 {
  596.                     word+=base;
  597.                     carry--;
  598.                 }
  599.                 a1[digit+j] = ((PlatWord)(word%base));
  600.             }
  601.         }
  602.         LISPASSERT(carry == 0);
  603.             
  604.         //D5:
  605.         aQuotient[j] = (typename T::ElementType)q;
  606.         //D7:
  607.         j--;
  608.         
  609.     }
  610.  
  611.     //D8:
  612.     a1.SetNrItems(n);
  613.     PlatDoubleWord carry;
  614.     BaseDivideInt(a1, d, base,carry);
  615.     aRemainder.CopyFrom(a1);
  616. }
  617.  
  618.  
  619.